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Systematic Study of Rogue Wave Probability Distributions in a 
Fourth-Order Nonlinear Schrodinger Equation 
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Abstract. Nonlinear instability and refraction by ocean currents are both important 
mechanisms that go beyond the Rayleigh approximation and may be responsible for the 
formation of freak waves. In this paper, we quantitatively study nonlinear effects on the 
evolution of surface gravity waves on the ocean, to explore systematically the effects of 
various input parameters on the probability of freak wave formation. The fourth-order 
current-modified nonlinear Schrodinger equation (CNLS 4 ) is employed to describe the 
wave evolution. By solving CNLS numerically, we are able to obtain quantitative pre- 
dictions for the wave height distribution as a function of key environmental conditions 
such as average steepness, angular spread, and frequency spread of the local sea state. 
Additionally, we explore the spatial dependence of the wave height distribution, asso- 
ciated with the buildup of nonlinear development. 



1. Introduction 

Oceanic rouge waves, or freak waves, are surface grav- 
ity waves whose wave heights are extremely large compared 
to the typical wave in a given sea state. Freak waves have 
been reported throughout maritime history, as they are a 
hazard to mariners, cargo ships, and even to large cruise 
liners. However, serious scientific investigation commenced 
more recently, due to the availability of wave state mea- 
surement methods including satellite remote sensing, and 
the introduction of stochastic process theory to ocean wave 
forecasting [Bates, 1952; Kinsman, 1965]. 

The scientific community has studied the topic experi- 
mentally (in water tanks) [Onorato et al, 2004] and [Ono- 
rato et al., 2006], observationally (with satellite imaging, 
for example) [Dankert et al, 2003] and [Schulz-Stellenfleth 
and Lehner, 2004], and theoretically [Kharif and Pelinovsky, 
2003]. A review of recent progress in freak wave forecasting 
can be found in an article by Dysthe et al. [2008] . Given the 
chaoticity of ocean wave dynamics, it is not possible to pre- 
dict individual freak wave events. Instead, we study wave 
behavior probabilistically. The simplest theory assumes a 
random superposition of a large number of monochromatic 
waves with different frequencies and propagating directions, 
as in the Longuet-Higgins random seas model [Longuet- 
Higgins, 1957]. By the central limit theorem, this model 
leads to a Rayleigh probability distribution of wave heights, 
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where 2H is the wave height, a 2 is the variance of the sur- 



face elevation, and V27tct is the mean wave height. For linear 
waves, the wave height 2H (crest to trough) equals twice the 
crest height due to a symmetry of the wave equation. This 
symmetry is broken for nonlinear waves, and throughout 
the paper, we focus on the wave height 2H , where the crest 
height is H to leading order. 

Conventionally, freak waves are defined as having wave 
height 2H > 8.8ct, or alternatively 2H > 2.2H S , where H s 
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is the significant wave height [Dean, 1990]. The significant 
wave height H s was traditionally defined as the average of 
the tallest one-third of waves in a wave train, which comes to 
H a « 4.004cr when the height distribution is Rayleigh. Since 
the Rayleigh distribution (1) serves as an excellent zeroth- 
order approximation for realistic wave height distributions 
in the ocean, H a = 4ct may be used today as another defini- 
tion of the significant wave height [Forristall, 1978]. Simi- 
larly, extreme freak waves are defined as having wave height 
2H > 12. Oct or 2H > 3H S . According to the simple Rayleigh 
distribution (1), the occurrence probabilities of a freak wave 
and an extreme freak wave are 6.3 x 10 -5 and 1.5 x 10 -8 , re- 
spectively. Compared with the observational data [Dankert 
et al., 2003], the stochastic Rayleigh model greatly under- 
estimates the probability of freak waves. Recent advances 
in technology have allowed multiple wave tank experiments 
and field observations to be conducted, confirming the need 
for a more realistic theory to explain the results [Forristall, 
2000; Onorato et al, 2004]. 

Nonlinear instability and refraction by ocean currents are 
both important mechanisms that go beyond the Rayleigh 
approximation and may be responsible for the formation of 
freak waves. 

The nonlinear instability, also known as the Benjamin- 
Feir instability in the one-dimensional case [Benjamin and 
Feir, 1967], describes how a monochromatic wave can be 
unstable under a class of small perturbations. In the two- 
dimensional case, however, the spread in the wave vector 
happens not only in magnitude but also in wave direc- 
tion. The nonlinear evolution involves both a nonlinear self- 
focusing effect at short times, which leads to an increase in 
the occurrence rate of extreme events, and also a disper- 
sion in wave vector on larger time scales, which leads to a 
decrease in the probability of extreme events as the initial 
spectrum becomes less unidirectional [Onorato et al., 2002]. 

We usually use the steepness to describe t he n onlinearity 
of a water wave. Steepness is defined as e — 2Hka/2, where 
2H is the mean wave height and fco is the peak wave number. 
The nonlinear evolution of a surface-gravity water wave is 
believed to be well described by the n onli near Schrodinger 
equation when the steepness is small, 2Hko <C 1, and the 
wave bandwidth is narrow, Afc|/fco ~ s <C 1. This is sup- 
ported by Trulsen and Dysthe's analysis of the famous New 
Year's wave event, which occurred at the Draupner platform 
in the North Sea in 1995 [Trulsen and Dysthe, 1997]. 

Clearly the exact evolution is always nonlinear to some 
extent, but the key is to introduce nonlinearities at the right 
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moment and in an insightful and computable way. Realis- 
tic fully nonlinear computations wave by wave over large 
areas are very challenging, but initial attempts have been 
made to simulate the ocean surface using the full Euler equa- 
tion both on large scales [Tanaka, 2001] and over smaller 
areas [Gibson and Taylor, 2005; Gibson and Swan, 2007]. 
However, in this paper, we rely on series expansions, re- 
sulting in a sequence of approximations. Pioneering work 
using a third-order equation in e (the original Nonlinear 
Schrodinger Equation (NLS)) was performed by Zakharov 
[1968], and then extended to fourth order by Dysthe [1979] 
(for water waves propagating in an ocean of infinite depth). 
In nature, the small steepness constraint is usually satisfied, 
since the steepness under storm conditions is normally less 
than 0.09 [Dysthe et al, 2008]. However, the narrow band- 
width constraint is often not complied with in naturally oc- 
curring conditions. To overcome this problem, a modified 
nonlinear Schrodinger equation (MNLS) has been obtained 
by Trulsen and Dysthe [1996] for broader bandwidth, only 
requiring jAfc|/fc ~ \fe, e«l. 

The effect of currents interacting with an incoming sea 
state has been studied by Heller et al. [2008] using ray 
dynamics and by Janssen and Herbers [2009] using Monte 
Carlo simulations. These studies suggest that the interac- 
tion between incoming wave and current may serve as a 
triggering mechanism for the formation of freak waves dur- 
ing nonlinear wave evolution. In the work by Heller et al. 
[2008], a freak index is defined in terms of the mean wave 
speed, mean current speed, and the angular spread of the 
incoming wave, and a quantitative relationship is predicted 
between this freak index and the occurrence probability of 
freak waves. However, in this work ray dynamics is used in 
place of the real wave equation for ocean waves, and nonlin- 
earity is not included in the model. The calculations done 
by Janssen and Herbers [2009] do include nonlinearity, and 
serve as a proof of principle for the importance of the in- 
teraction between nonlinearity and refraction by currents. 
Still, results are obtained only for specific values of the in- 
put parameters. Thus, a systematic quantitative analysis 
of the freak wave probability as a function of the strength 
of nonlinear and refractive effects has not previously been 
undertaken. 

Here, we use the current-modified nonlinear Schrodinger 
equation (CNLS) developed by Stocker and Peregrine [1999] 
to describe the evolution of the wave envelope, with the aim 
of integrating both nonlinear effects and wave-current inter- 
action in one model. This approach allows for a full quan- 
titative analysis of both these mechanisms of freak wave 
formation. In the present paper we focus on numerically 
exploring the dependence of the wave height distribution on 
nonlinear effects only, by setting the currents to zero. We 
will see that interesting behavior is obtained already in the 
current-free regime. The joint dependence on steepness and 
current strength for nonzero random currents will be pre- 
sented in a future publication. 

This paper is organized as follows: In section 2 we re- 
view the CNLS equations and present the basic computa- 
tional setup. In section 3 we discuss the expected form of 
the wave height distribution in the presence of deviations 
from Rayleigh statistics due to nonlinearity and currents. 
In section 4 we systematically study the dependence of the 
wave height distribution on the parameters describing the 
incoming sea state, allowing the probability of freak wave 
formation to be calculated as a function of these param- 
eters. Finally, in section 5 we summarize our results and 
discuss the outlook for the future. 

2. Model 



The current-modified nonlinear Schrodinger equation 
equation (CNLS) extends the NLS to include a random sur- 
face current varying around a mean value uo- The time- 
independent current field can be expressed as: 

u(r) =u + 8u(r) , (2) 

where the random current velocity fluctuations 5u are as- 
sumed to be 0(e 2 ), and slowly varying on the scale of a 
wavelength. For surface gravity waves in deep water (i.e., 
when the water depth is much larger than the wavelength), 
the dispersion relation is given by: 

oj(k,f) = ^Jg~$\+k-u{r), (3) 

where uj is the angular frequency of the incoming wave whose 
wavelengt h is g iven by A = 2ix/k, and the wave velocity is 
VfcW = \J g/4k k + u. For convenience, we take the peak 
wave vector of the incoming wave to be in the positive x 
direction. The wave vector of each wave component is then 
expressed as: 

k = k x + (dk x ,dky) , (4) 

where fco > 0. 

To simplify the equations and simulations, we work in 
th e fram e of reference moving with the peak wave velocity 
\J g/4:ko x + uo, and study the evolution of the wave enve- 
lope instead of the wave function itself. The complex wave 
envelope a(r, t) is defined by 

C(f, t) = Re a(f, t y k «*-*V^t > (5) 

where £(r, t) is the surface elevation. For a narrow-banded 
spectrum with small angular spread (i.e., |A/c|/feo ~ y/e <C 
1), the wave envelope varies on a scale long compared to 
the wavelength. In that regime, the probability distribution 
of the wave envelope a(r, t) closely approximates the crest 
height distribution. 

Finally, it is co nvenie nt to employ dimensionless variables 
A — koa, b"U — yJko/g8u, and 

(X, Y, Z, T) = (k x ~ ^M,k y, k z, yfgT t) . (6) 

In these variables, the third-order CNLS equation is [Stocker 
and Peregrine, 1999]: 

iA T -lAxx + ]A Y Y-l;A\A\ 2 -ASUx = 0, (7) 

o 4 2 

where subscripts denote partial derivatives. For larger steep- 
ness, i.e., larger A, the third-order equation does not suffi- 
ciently account for nonlinear effects, and we need the fourth- 
order CNLS (CNLS 4 ) equation, also given by Stocker and 
Peregrine [1999]: 

iB T - l(B X x - 2B YY ) - \b\B\ 2 - B$ cX 
= ^.{Bxxx - 6B YYX ) + ^B(BB* X - 6B*B X ) 
+$xB + i(^cxr - $cz)B - iVn $ c • V H B . (8) 

On the left hand side of (8) are the third-order terms (for 
comparison with (7)), while the right hand side contains all 
terms appearing at fourth order. Here Vh_= (d/dX,d/dY) 
is the gradient in the horizontal plane, and $, $ c , and B rep- 
resent the mean flow, surface current, and oscillatory parts, 
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respectively, of the velocity potential 




Re (Be kaz+i,fi + B 2 e 2( *°* +<rt )] ,(9) 



where ip = kox — y/gkot — X — T/2 is the phase. The surface 
elevation, which is the quantity of interest for our purposes, 
is similarly expanded as 

C^C + Cc + fco" 1 [Re (Ae* + A 2 e 2lv + A 3 e 3 ^)] ,(10) 

where the expansion coefficients may be obtained from the 
velocity potential as 



was conducted on a 512 by 1024 grid, corresponding to a 10 
km by 20 km ocean area or 64 by 128 wavelengths. 

The initial incoming wave is a random linear superpo- 
sition of A/ - monochromatic plane waves with different fre- 
quencies and propagating directions: 
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Each wave vector fc, in (12) can be expressed as 
k = (ko + k') ■ (cos 6 x + sin y) , 



(12) 



(13) 



.4 
A 2 
A 3 



= iB 



2k 



B x + 



2B. 



yy) + g B \ B \ 



-is 2 + ±BB a 
2 k 

— -B 3 . 

8 



(11) 



We have performed calculations using both third-order 
and fourth-order equations, but all results shown below are 
obtained using CNLS 4 . A split -operator Fourier transform 
method developed by Weidman and Herbst [1986] was em- 
ployed to solve the PDE numerically. The evolution equa- 
tion is separated into two parts: a free evolution part and 
a part containing nonlinear and current terms. In each 
time step, the wave envelope is transformed into momen- 
tum space where the free evolution operator is applied, and 
then transformed back into position space and acted on by 
the nonlinear and current operator. In the case of CNLS , 
the free evolution part is: 



-(Bxx 



2Byy) + —{Bxxx 



GByyx) 



and the nonlinear and current part is: 



1 , ,■ 
~ 2 B\B\ 



B$ cX + $xB + -B{BB X ~ QB'Bx) 



+i{^cXT - §cz)B - iV H $c ■ VhB . 

Without loss of generality, we choose the mean wave- 
length of the incoming wave to be 156 m, corresponding 
to a mean wave speed of 7.8 m/s. Considering both the 
computation time and statistical accuracy, the calculation 
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Current-Modified Nonlinear Schriidin 
Equation (CNLS) 



Initialization aiea Wave development aiea Absorbing area 

Figure 1. A sketch showing the calculation setup. The 
10 km by 20 km area is discretized using a 512 by 1024 
grid. The incoming wave is prepared at the left side of 
the grid, developing along the x axis from left to right, 
and is eventually absorbed in the right boundary. 



where the random wave number variation k' follows a Gaus- 
sian distribution of half height width Afc, and the angle 
6 from the mean propagation direction is uniformly dis- 
tributed in the interval [— /3A0, V3A0]. (To leading or- 
der, the results depend only on the standard deviation, A8, 
and not on the higher moments of the angular distribution.) 
Each coefficient a is drawn from a normal distribution, 
and the final sum is normalized to have the desired mean 
wave height. We have confirmed that our results are Af- 
independent for sufficiently large AT; specifically, AT = 2 • 10 5 
is used for obtaining the data shown in this paper. 

The basic setup of the calculation is shown in Figure 1. 
The incoming wave is prepared as a random superposition 
of plane waves at the left side of the grid, developing along 
the x axis from left to the right, and is eventually absorbed 
at the right side boundary. The split-operator Fourier trans- 
form method used in the calculation automatically imposes 
periodic boundary conditions on the system. Absorption 
of the outgoing wave on the right is required to prevent it 
from re-entering the grid on the left side. However, to avoid 
reflection or other boundary effects, we allow the outgoing 
wave on the right to decay gradually using a tanh() multi- 
plicative factor, and similarly another tanh() multiplicative 
factor allows the incoming wave to build up gradually in the 
initialization area on the left side of the grid, as indicated 
by the dotted lines in Figure 1. The boundary conditions 
in the transverse (y) direction remain periodic. By varying 
the system size in both x and y, we have confirmed that 
for the parameters used below, boundary conditions have a 
negligible effect on the wave height statistics. 

3. Form of the Wave Height Distribution 

In the Longuet-Higgins random seas model [Longuet- 
Higgins, 1957], the sea state is given by a random super- 
position of many plane waves with differing directions and 
frequencies, and by the central limit theorem, the surface 
elevation function £(r) is distributed as a complex Gaus- 
sian random variable with standard deviation a. Further- 
more, for a narrow-banded spectrum with small angular 
spread (|(5fc| <C |fc|), the crest height H is equal to the wave 
function amplitude |£|, and the probability distribution of 
wave height 2H is given by the Rayleigh distribution of (1). 
The corresponding cumulative probability of encountering a 
wave of height 2H or larger is 



-Paayleigh(2-ff) = 6 



-(2H) J /(8(T 2 ) 



(14) 



We now consider the effect of nonlinearity and currents, 
where the wave envelope a(f, i) defined by (5) is evolving 
in accordance with the equations presented in the previous 
section. Since we still assume a narrow spectrum and a 
small angular spread, the envelope evolves slowly in space 
and time, on the scale of the mean wavelength and mean 
wave period, respectively. In analogy with the situation 
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for linear ray dynamics in the presence of random cur- 
rents [Yvng et al., 2011], in the neighborhood of any space- 
time point (r.t), we have a wave intensity proportional to 
H 2 — \a(r,i)\ = |C(f, t)\ 2 . Thus, in contrast with the orig- 
inal Longuet-Higgins model, the locally-averaged wave in- 
tensity is not uniform over all space and time but instead 
exhibits "hot spots" and "cold spots" in space-time associ- 
ated with energy focusing and defocusing. Hence the wave 
height distribution is given locally in space and time by a 
Rayleigh distribution around the local mean height (corre- 
sponding to a locally random superposition of plane waves), 
while the local mean height itself varies slowly on the scale 
of the mean wavelength and mean period. 

At each point in space-time, the central limit theorem 
and thus the Rayleigh distribution still apply, and we have 



P^ t) (2H) = e 



~(2Hf/(8v 2 I(r,t)) 



where I(r, t) is the local energy density, normalized by the 
average energy density of the incoming sea state (so that 
1 = 1), and a is the variance of the surface elevation in the 
incoming sea state, before the nonlinear evolution begins to 
act. 

Now averaging over space and over time, we obtain a total 
cumulative wave height distribution 



P taM (2H)= / dlg(l) 



-(2H) 2 /(8<r 2 /) 



(16) 



In (16), the full cumulative distribution of wave heights for 
a given sea state has been expressed as a convolution of two 
factors: (i) the local density distribution g(I), which at least 
for small deviations from Rayleigh is conjectured to be \ 2 
distributed in analogy to the ray dynamics approximation, 
and (ii) the universal Longuet-Higgins distribution of wave 
heights for a given local density. Similar decompositions of 
chaotic wave function statistics into non-universal and uni- 
versal components have found broad applicability in quan- 
tum chaos, including for example in the theory of scars [Ka- 
plan, 1999; Smith and Kaplan, 2009; Backer and Schubert, 
2002]. In the context of rogue waves, a similar approach was 
adopted by Regev et al. to study wave statistics in a one- 
dimensional inhomogeneous sea, where the inhomogeneity 
arises from the interaction of an initially homogeneous sea 
with a (deterministic) long swell [Regev et al., 2008]. 

Taking the local mean intensity to be \ 2 distributed with 
N degrees of freedom, 



in the same dimensionless units. 

Notably, the N parameter describing deviations from 
Rayleigh statistics may be directly _related to the excess 
kurtosis of the sea state, 72 = V 4 /(v 2 ) 2 ~ 3, which mea- 
sures deviations from Gaussian statistics for the surface 
elevation r\ [Onorato et al, 2002; Janssen and Herbers, 
2009]. For a narrow-banded spectrum, we have 72 = 

(3/2)(2J/) 4 /(2ff)2 2 - 3, and from (18) we obtain 



72 = 



_6_ 

N 



(21) 



4. Numerical Results 



(15) 4.1. Wave Height Probability Distributions 



In our simulations, we choose the mean wave number to 
be ko = 40 km -1 (a typical value in a normal sea state), 
corresponding to a group velocity vo = 7.8 m/s, wavelength 
Ao = 156 m, and period To = 10 seconds. Note that there 
is no loss of generality in this choice, as it merely sets the 
fundamental spatial and time scales for the wave dynam- 
ics. Each run simulates wave evolution for 4 ■ 10 5 seconds, 
corresponding to 4 • 10 4 periods. In the setup pictured in 
Figure 1, the wave front takes approximately t = 2000 sec- 
onds to travel from the initialization area to the absorb- 
ing area taking the straight line path, and the system fully 
equilibrates at around t = 7000 seconds. The time interval 
7000 s < t < 400000 s is used to obtain wave height statistics 
in the wave evolution region defined by 4 km < x < 16 km, 
which indicated by the shaded area in Figure 1. Al l w ave 
heights are normalized by the mean wave height 2H in 
a given simulation run, obtained using the time interval 
7000 s < t < 20000 s. In the end, we calculate the prob- 
ability distribution of wave heights in each run, and by re- 
peating this process for different input parameters we obtain 
the dependence of the wave height distribution on the sea 
parameters. 

At present, we set the current to zero, and focus on the 
nonlinear effects. In this case, we find that the wave height 
distribution for given sea parameters is reliably obtained 
from a single run given a sufficiently long run time, as dis- 
cussed above. The two important properties of the incom- 
ing wave are the steepness and the wave vector spread. The 
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(17) 



and convolving the \ 2 distribution of the mean intensity 
with the Rayleigh distribution around the mean intensity, 
we can obtain as in the linear case [Hohmann et al., 2010] a 
K-distribution for the total distribution of wave heights: 



Ptotal(H) = 2- 



NH/2a) ' 



T(N/2) 



-K 



N/2 



NHa) 



(18) 



where K„(y) is a modified Bessel function. 

Defining the dimensionless variable x = 2H/H S ~ 
2H/{Aa), so that a rogue wave is given by x = 2.2 and 
an extreme rogue wave by x = 3.0, we find the probability 
of a wave height exceeding x significant wave heights: 



(VNx)' 



-K N/2 (2VNx) 



T(N/2) 

to be compared with the random seas prediction 

D I \ —2a; 2 

"Rayleigh l,^ J — C 



(19) 



(20) 




Wave Heighi/SWH 



Figure 2. The distribution of wave heights, in units of 
the significant wave height, is calculated for four nonzero 
values of the steepness e (upper three dashed or dotted 
curves), and compared with the random seas model of 
(20) (lowest solid curve). In each case, the solid curve is 
a best fit to the K-distribution of (19). Here the we fix 
the angular spread AO = 2.6° and wave number spread 
Afc/fco = 0.1 of the incoming sea. 
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steepness e = ko2H /2 controls the relative importance of 
the nonlinear term i-B|B| 2 in (8) compared with the linear 
term s; the steepness is adjusted by varying the mean height 
2H of the incoming sea, since the mean wave number fco is 
fixed throughout. 

Two parameters are needed to describe the wave vector 
spread: the variation in magnitude (wave number) is charac- 
terized by the ratio Afc/fco (or equivalently by the frequency 
spread Aui/uq ~ (l/2)Afc/fco), and the angular spread is 
given by A9. The objective is to find the response of the 
freak wave occurrence probability to these three input pa- 
rameters. 

First, we keep the wave vector spread constant, and vary 
the steepness in the range 0.016 < e < 0.07. The steep- 
ness as defined above is normally less than 0.09 under storm 
conditions [Dysthe et al, 2008]. Since we are focused on 
the normal sea state, ignoring for example tsunamis (large 
tidal waves caused by underwater earthquake), we are not 
interested in the behavior at anomalously large values of the 
steepness. The lower limit is determined by two consider- 
ations. Physically, when the steepness becomes very small, 
the distribution of wave heights approaches the Rayleigh 
limit, where the wave height distribution is already well un- 
derstood, and the number of extreme waves is very small. 
Numerically, at small values of the steepness, the finite simu- 
lation time makes it difficult to quantify deviations from the 
baseline Rayleigh distribution, due to the paucity of events 
in the tail of the wave height distribution. 

Typical results are represented by dashed or dotted curves 
in Figure 2, where we fix Ak/k = 0.1 and AO = 2.6° (the 
values of A6 required to see very strong effects from non- 
linear focusing are typically smaller than those needed to 
observe significant deviations from Rayleigh by linear scat- 
tering [ Ying et al., 2011]). The cumulative probability distri- 
bution of the wave height 2H , in units of the significant wave 
height H s , is shown for three nonzero values of the wave 
steepness e. From bottom to top, the three thick dashed or 
dotted lines show results for steepness e = 0.019, 0.032, and 
0.042. As expected, the Rayleigh probability distribution 
of (20) is recovered in the limit e — > 0, and ever stronger 
enhancement in the tail is observed as the steepness of the 
incoming sea increases. The occurrence probability of ex- 
treme rogue waves, 2H/H S — 3.0, is enhanced by one to 
three orders of magnitude for the parameters shown. 

The shape of the probability distribution and the typ- 
ical enhancements in the probability compared with the 
Rayleigh baseline are consistent with the results of previous 
numerical simulations. For example, Onorato et al. [2006] 
have obtained an enhancement of up to one order of magni- 
tude in the probability of freak wave formation (2H/Hs > 
2.2) and up to 2.5 orders of magnitude in the probability of 
extreme freak wave formation (2H/H S > 3), for e = 0.08 
(k H a /2 = 0.125), Afc/fco = 0.18 (Alj/lu = 0.09), and 
AO = 0. 

In Figure 2, each data set is fit to the K-distribution of 
(19), arising from the local Rayleigh approximation. We see 
that the fits, indicated by solid lines, perform adequately 
for probabilities down to 10~ 6 , where statistical noise be- 
gins to dominate. The best-fit values are N = 53.2, 12.4, 
and 8.2, with smaller steepness corresponding to larger N. 
The Rayleigh distribution (N — > oo) is also shown here, 
corresponding to e = 0. 

In particular, we clearly observe the crossover between 
the Gaussian behavior (20) at small to moderate heights 
and asymptotic exponential behavior at large heights. How- 
ever, systematic deviations do exist, which are especially 
visible at larger values of e, corresponding to smaller val- 
ues of the N parameter. These systematic deviations are in 
large part due to the fact that the true wave height distri- 
bution for any given set of input parameters exhibits spatial 
dependence, evolving from the original Rayleigh distribution 



imposed by incoming boundary conditions to the broader re- 
distribution, and then gradually back to a Rayleigh distribu- 
tion as the wave energy is transferred to longer wavelengths 
and the steepness decreases [Janssen and Herbers, 2009]. 
An example of this spatial dependence appears below in Fig- 
ure 3. Thus, a more accurate model for the total wave height 
distribution consists of a sum of several K-distributions, or 
equivalently the tail of the full distribution may be modeled 
by a K-distribution multiplied by a prefactor C < 1. Nev- 
ertheless, as seen in Figure 2, (19) correctly describes wave 
height probabilities at the ±20% level of accuracy, allows 
for an extremely simple one-parameter characterization of 
the wave height distribution, and facilitates easy compari- 
son between the effects of linear and nonlinear focusing. 

4.2. Effect of Several Sea Parameters on the Wave 
Height Distribution 

Based on Figure 2, we can clearly qualitatively conclude 
that larger steepness will lead to a higher freak wave occur- 
rence probability. Now by fitting the distribution to (19), 
in Figure 4, we examine quantitatively how increasing the 
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Figure 3. The fourth moment (2H) 4 of the wave height 
distribution, in units where 2H = 1, is shown as a func- 
tion of evolution distance for different values of steepness, 
starting in each case from a Longuet-Higgins random 
sea with mean wave speed v — 7.81 m/s, initial angular 
spread A9 — 1°, and wave number spread Ak/kg = 0.1. 
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Figure 4. The best-fit N value in (19) describing the 
wave height probability distribution is shown as a func- 
tion of the steepness e for several values of the initial 
angular spread A6. The initial wave number spread is 
fixed at Ak/ko = 0.1, as in Figure 2. The line shows the 
best-fit scaling N ~ e" 2 ' 9 for AO = 2.6°. 
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steepness will affect the N value. In the figure, each point 
corresponds to one particular sea state, and points of a given 
shape correspond to the same initial angular spread but dif- 
ferent values of the initial steepness. The horizontal axis 
indicates the steepness, and the vertical axis shows the A^ 
value. From bottom to top, the four curves correspond to 
angular spread AO = 1°, 2.6°, 3.6°, and 5.2°. A power 
law is observed on the log-log scale for each angular spread. 
Taking AO — 2.6° (indicated by squares in Figure 4) as an 
example, we obtain the scaling 



N . 



(22) 



where c ~ —3. 

The AO = 3.6° and A9 — 5.2° curves lie very close to 
each other, indeed, at larger values of the steepness (not 
shown), saturation occurs, indicating that the occurrence 
rate of freak waves does not fall below a certain level when 
the steepness is sufficiently high, independent of the angle 
AO. However, we also note that the CNLS equations are 
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Figure 5. The horizontal axis indicates the initial an- 
gular spread AO, and the vertical axis shows the best-fit 
N value, which determines the wave height probability 
distribution. Here the steepness is fixed at e — 0.032, 
and results are shown for four values of the wave number 
spread Afc/fco. The best-fit line shows the scaling of N 
with AO in one example, specifically N ~ (Ad) 1 ' 04 for 
Ak/ko = 0.15. 




Figure 6. The horizontal axis indicates the wave num- 
ber spread Ak/ko, and the vertical axis shows the best-fit 
N value, which determines the wave height probability 
distribution. Again, the steepness is fixed at e — 0.032. 
The best-fit line shows the scaling of N with Ak/ko, m 
this example, N ~ (Afc/fco) 1 ' 15 for AO = 2.6°. 



based on a small steepness approximation, and thus the be- 
havior obtained at very high values of the steepness may not 
be reliable. In the more realistic low-to-moderate steepness 
regime shown in the figure, N grows with increasing initial 
angular spread, i.e., the occurrence rate of freak waves falls. 

Now that we have an understanding of how A^ changes 
with steepness, the next step is to fix the steepness and vary 
the initial angular and wave number spread. At the end, we 
will be able to bring these results together, to obtain the 
dependence of the freak wave occurrence rate on all three 
key environmental parameters. 

Figure 5 shows the dependence of A^ on the initial angular 
spread for several different values of the frequency spread, 
while the steepness is fixed at 0.032. In this figure, points 
of a given shape represent sea states with the same wave 
number spread. The initial angular spread is varied on the 
horizontal axis. 

From bottom to top, the four data sets represent wave 
number spreads Ak/ko = 0.2, 0.15, 0.1, and 0.08, and the 
dependence on initial angular spread is well represented by 
N ~ (A0) a , where a = 0.96, 1.04, 1.07, and 1.03, for the 
four data sets, respectively. Thus, N scales as the first power 
of the angular spread at this value of the steepness, indepen- 
dent of the wave number spread. 

Similarly, we can choose several typical values of the ini- 
tial angular spread, and obtain the dependence of N on the 
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Figure 7. The best-fit N value in (19) describing the 
wave height probability distribution is shown as a func- 
tion of wave vector spread (keeping Ak/ko = 7 • A0/n 
while varying Afc/fco for several values of the steepness 
e). The line shows the best-fit scaling with power 2.32 
for e = 0.032. 
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Figure 8. The angular spread AO is varied for several 
values of the steepness e, keeping the wave number spread 
Afc/fco = 0.1 fixed. From bottom to top, the four data 
sets correspond to e = 0.07, 0.042, 0.032, and 0.026. 
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initial wave number spread, as shown in Figure 6. Again, 
the steepness is fixed at 0.032, and the other physical condi- 
tions are the same as in the previous analysis. All sea states 
having the same initial angular spread are represented by 
the same symbol. From bottom to top, the angular spread 
is A9 — 1°, 2.6°, 3.6°, and 5.2°. All the initial sea states in 
this figure have steepness e = 0.032. The power law is con- 
sistent for different initial angular spread, showing that N 
scales approximately as the first power of the wave number 
(or frequency) dispersion. 

We notice deviations from the power-law behavior near 
the left edge of Figures 5 and 6, corresponding to very 
small angular spread A9 w 1° and very small wave num- 
ber spread Ak/ko ~ 0.08 (Aw/wq w 0.04), respectively. At 
even smaller values of A6 and Ak (not shown), the CNLS 
equations we rely on, which are based on a perturbative ex- 
pansion in powers of the surface elevation, begin to break 
down. This breakdown is evidenced, for example by large 
discrepancies appearing between the wave height distribu- 
tions obtained using the fourth-order CNLS (8) and third- 
order CNLS (7), which make the CNLS expansion untrust- 
worthy. Furthermore, we recall that the functional form of 
the K-distribution (18) is obtained assuming modest fluctu- 
ations of the local intensity around the mean intensity, i.e., 
this functional form makes sense only for N > 1, where 
the wave height distribution may be considered to be a 
modified Rayleigh distribution. We note in particular that 
the Benjamin-Feir modulational instability is kn own t o be 
present in the NLS equations for (Afc/fco)/s < ^32/n [Al- 
ber, 1978], and for sufficiently small values of {Ak/ko)/e, 
this instability will dominate the wave evolution so that ex- 
treme waves become commonplace rather than rare events. 
This regime is outside the range of validity of the present 
analysis. Nevertheless, as we will see below, the moderate 
parameter values for which our analysis is applicable, which 
are also the parameters most likely to occur in nature, can 
given rise to enhancement factors as large as 10 3 in the like- 
lihood of occurrence of extreme freak waves (Table 2). 

4.3. A Unified Model 

In the previous subsection, we examined separately the 
dependence of the wave height distribution on steepness, 
angular spread, and wave number spread, keeping the other 
two variables fixed. We would now like to understand the 
combined dependence on these three parameters describing 
the initial sea state. 

Since both angular and wave number spread control the 
area in wave vector space in which wave modes are per- 
mitted, we examine a scenario in which the angular spread 
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Figure 9. The initial wave number spread Ak/ko is 
varied for several values of the steepness e, keeping the 
initial angular spread A0 — 5.2° fixed. From bottom to 
top, the four data sets correspond to e = 0.07, 0.042, 
0.032, and 0.026. 



and wave number spread change at a constant ratio for a 
given value of the steepness. We arbitrary choose the ratio 
as Ak/ko = 7 • A9/ir. The results are shown in Figure 7, 
again plotted on a log-log scale, with the frequency spread 
Ak/ko on the horizontal axis. Each symbol represents all sea 
states with a given value of the steepness. The line shows 
the best-fit power law N ~ (Afc/fc ) 2 ' 32 for e = 0.032. From 
bottom to top, the other three curves correspond to steep- 
ness e — 0.051, 0.042, and 0.026, corresponding to a scaling 
with power 2.3, 2.0, and 2.5, respectively. 

First, we focus on the case of steepness e = 0.032, where 
we have a clear understanding of how N depends on the an- 
gular and frequency spread. We notice that the power law 
exponent 2.32 appearing in Figure 7 when Ak/ko and A6 
are varied together is close to the sum of the exponent ob- 
tained when Ak/ko and A9 are varied separately. Thus, we 
propose a model to describe the dependence of N on angular 
and frequency spread: 

»-»(£)'(")'■ 

where B, a, and ft may depend on the steepness. Now, if we 
fix the ratio between the angular and wave number spread, 
*M=A.£* (23) yields 




Comparing (23) and (24) with the data shown in Figure 5 
and Figure 7, and setting A — 1/7, we find the results are 
consistent and give B — 5279, a = 1, and ft — 1. Thus, we 
obtain the following predictive model (when e = 0.032): 

"-■»■(£)(¥)■ < a » 

The next step is to find out how B, a, and /3 in (23) may 
evolve when the steepness changes. For several typical val- 
ues of the steepness e, we vary the initial angular spread A9, 
while keeping the frequency spread fixed at a typical value, 
Ak/ko = 0.1. The results are shown in Figure 8. Here each 
symbol represents sea states with one value of the steepness: 



Table 1. Relation Between the Steepness e and the Coeffi- 
cient B in (26) 



Steepness e 


B 


0.026 


17700 


0.032 


5000 


0.042 


1600 


0.052 


1100 



Table 2. Enhancement in the Probability of Rogue Wave 
Formation a 



N 


E(2.2) 


£(3.0) 


72 


2 


1.1 ■ 10* 


5.2 • 10 4 


3 


5 


37 


7.3 ■ 10 3 


1.2 


10 


16 


1.3 • 10 3 


0.6 


20 


6.8 


2.2 ■ 10 2 


0.3 


50 


2.9 


27 


0.12 


100 


1.8 


7.8 


0.06 



a The enhancement in the probability of rogue wave forma- 
tion (wave height 2H = 2.2 H s ) as well as the enhancement of 
the probability of extreme rogue wave formation (wave height 
2H = 3.0 H s ) and the excess kurtosis 72 (defined in (21)) are 
calculated for several values of the N parameter. Here E(x) = 
ftotal(aO/PB.ayleigl»(aO- 
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from bottom to top, the data sets correspond to steepness 
e = 0.051, 0.042, 0.032, and 0.026. Fitting each data set 
to a power law relating N and A9, we find that the power 
(a in (23)) is 0.93, 1, 1.15, and 1.19 for these four values 
of e, respectively. Thus, a is negatively correlated to the 
steepness, i.e., the dependence on angular spread becomes 
slightly weaker as the steepness increases (consistent with 
the saturation effect noted earlier). Nevertheless, the effect 
is quite small for moderate steepness values, so to a good 
approximation we may treat a as a constant near unity. 

Similarly, in Figure 9 we fix the initial angular spread 
fixed at a typical value A9 = 5.2°, and examine the N de- 
pendence on the wave number spread Afc/fco for different 
values of the steepness e. Each wave height probability dis- 
tribution is fitted to obtain the N parameter, and a power 
law relating N and Afc/fco (the /3 exponent in (23)) is ob- 
tained for each value of e. From bottom to top, the data sets 
correspond to steepness e = 0.052, 0.042, 0.032, and 0.026, 
and yield best fit exponents ft — 1.07, 0.98, 1.03, and 1.39, 
respectively. Since the variation in ft is non-systematic and 
weak, we conclude that /3, like a, may be approximated as 
an £-independent constant near unity. 

Thus, we conclude that in (23), both a and ft are inde- 
pendent of steepness s, and only B is a function of steepness: 



U A J 



(26) 



In Table 1, we present the value of B corresponding to val- 
ues of the steepness e. From Figure 4 and (22), we conclude 
that that N ~ e -3 , If we use this scaling to fit the data in 
Table 1, we obtain B = 0.203e -3 , and thus the final model 
is described by: 



™-'(t)(v) 



(27) 



As noted previously, the scaling breaks down for very small 
Afc/fco or A9, or very large e, for example, when the two- 
dimensional Benjamin-Feir index [Mori et al, 2011] is large. 

We now examine the implications for the occurrence prob- 
ability of a freak wave. According to the conventional def- 
initi on o f a frea k w ave, the wave height for a freak wave is 
3.51 2H , where 2H is the mean wave height. By integrat- 
ing the wave height distribution over heights larger than 
3.51 2H, we obtain the total freak wave occurrence proba- 
bility. For example, for e = 0.032, A9 = 5.2°, Ak/k = 0.1, 
we have N = 21.4. Comparing the result with a Gaus- 
sian theory, the freak wave occurrence probability is 6 times 
higher. Table 2 below shows the enhancement rate of the 
probability of freak wave occurrence in the nonlinear the- 
ory as compared with the simple Gaussian theory. We note 
that even at N values between 50 and 100, corresponding 
to the upper range of values in Figures 4 through 9, the oc- 
currence of extreme rogue waves is enhanced by an order of 
magnitude. Exponentially larger enhancement is predicted 
for parameters associated with smaller values of N. 

4.4. Spatial Dependence and the Kurtosis 

The distribution data in Figures 2 to 9 is taken over the 
entire area of wave development thus does not give informa- 
tion about the spatial build up of rogue waves. To study 
the spatial distribution of the wave heights, we can divide 
the wave development area into slices of thickness around 
0.3 km (approximately two wavelengths) in the wave prop- 
agation direction, and calculate the fourth moment of the 
wave height distribution, defined as (2H) 4 , in units where 
2H = 1, for each slice. For zero steepness, the waves are 
linear, and the wave height distribution is Rayleigh, giving 
a value of 3.242 for the fourth moment of the wave heights. 
Figure 3 shows the fourth moment of the wave height as 



a function of position for different values of the steepness, 
where the angular spread and wave number variation are 
fixed at A9 = 1° and Ak/k = 0.1. 

We can see from Figure 3 that the fourth moment of 
the wave height grows initially over a distance scale of sev- 
eral wavelengths, and then returns to the value given by 
a Rayleigh distribution. The result agrees with the Monte 
Carlo study conducted by Janssen and Berbers [2009]. 

In the following, we take a closer look at the wave height 
distribution in each slice. In this case we choose an exam- 
ple corresponding to the solid red curve in Figure 3, with 
AO = 1°, Afc/fco = 0.1, and steepness e = 0.054. To re- 
duce statistical noise, we divide the wave development area 
into larger slices of thickness around 1.3 km (approximately 
8 wavelengths). The cumulative probability distribution of 
the wave height 2H, in units of the significant wave height 
H 3 , is shown in Figure 10 for three different areas. The 
uppermost dot-dashed curve shows the result for the slice 



0.001 




E 
u 



2 3 4 

Wave Height/SWH 

Figure 10. The wave height distribution, in units of 
the significant wave height, is calculated for three dif- 
ferent spatial regions (dashed or dot-dashed lines), and 
compared with the Rayleigh random seas model of (20) 
(lowest solid curve). In each case, the solid curve is a 
best fit to a single K-distribution of (19). Here the we 
fix the angular spread A 6* = 1°, steepness e = 0.054, and 
wave number spread Afc/fco = 0.1 of the incoming sea. 
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Figure 11. The excess kurtosis 72 = 

(3/2)(2Jf) 4 /(2#) 22 - 3 is plotted versus distance 
along the direction of wave development. The upper 
data set (squares) is obtained directly from simulation 
results, and the lower data set(circles) is calculated using 
(21) from a K-distribution fitted to the wave height 
distribution in each slice. The parameters are the same 
as in Figure 10. 
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centered at the 7 km position; from Figure 3 we see that 
this is the peak region for wave height fluctuations. The 
short-dashed line shows the result for the slice centered at 
16 km, which is near the end of the wave development area, 
where wave statistics are gradually reverting to the Rayleigh 
limit. Finally, the long-dashed line in the middle is the wave 
height distribution for the whole wave development area, as 
indicated in Figure 1. 

In Figure 10, each data set is fit to the K-distribution 
of (19), arising from the local Rayleigh approximation. As 
noted in section 4.1, when wave height data are collected 
over a large spatial field that includes some areas of very 
strong deviations from Rayleigh statistics and other areas 
where such deviations have not yet had an opportunity to 
develop, the full distribution may not be well approximated 
by a single K-distribution. That is exactly what we observe 
in Figure 10: the wave height distribution for a single slice 
is better approximated by the K-distribution than that for 
the whole area. 

We can also observe that the wave height distribution in a 
slice in the ending area is better approximated by a single K- 
distribution than the distribution in the peak area. This is 
due to the fact that the distribution is evolving much faster 
in the peak area. Expanding to second order in slice thick- 
ness Ax for large N, the deviation from a K-distribution due 
to finite slice thickness scales as (Ax/N) 2 ■ (dN/dx) 2 , where 
dN/dx is the rate at which N evolves with position in the 
middle of that slice. This is consistent with the observation 
in Figure 3: the size of each slice is the same, but the slope 
dN/dx is larger and N is smaller in the peak area. There- 
fore, the deviation of the wave height distribution from a 
K-distribution due to finite slice thickness is larger in the 
peak area than in the ending area. 

In order to further investigate the deviation of the exact 
wave height distribution from a K-distribution, we compare 
for each slice the excess kurtosis 72 calculated directly from 
the simulation with the excess kurtosis obtained by fitting 
the wave height distribution to a K-distribution and then 
using (21). Of course, the answers would agree exactly if 
the distribution within each slice were perfectly described 
by a single K-distribution. The results are presented in Fig- 
ure 11, which includes data from 8 slices, ranging from the 
peak area (7 km) to the ending area (16 km). The initial 
nonlinear focusing area between 4 km and 7 km is omitted 
because the N value is varying too dramatically within a 
slice, making a fit of the wave height distribution to a single 
K-distribution meaningless. In the figure, we note that from 
right to left the difference between the exact kurtosis calcu- 
lation and the analytic formula of (21) increases systemat- 
ically, but decreases in the leftmost slice, centered on the 
peak. If we compare with the solid curve in Figure 3, we see 
that \dN/dx\ increases systematically from right to left and 
then decreases in the peak area. That is consistent with the 
results in Figure 11, given the scaling ~ (Ax/N) 2 (dN/dx) 2 
of the error due to finite slice thickness Aa;. 

5. Conclusions and Outlook 

In this paper, we study the dependence of the wave height 
distribution on the steepness, initial angular spread, and fre- 
quency of the incoming sea. We use the CNLS 4 nonlinear 
wave equation to simulate the ocean wave development, and 
obtain quantitative predictions for the wave height distribu- 
tion as a function of steepness, and angular and frequency 
spread. We then fit the wave height probability distribu- 
tion to a K-distribution, governed by a single parameter 
N that quantifies deviations from the Rayleigh distribution 
predicted by a simple random seas model. Furthermore, we 
show that N may be related to the excess kurtosis. We find 



simple power laws for the dependence of the N parameter 
on the environment parameters as long as the sea conditions 
are not extreme, and as a result, we propose a simple model 
to predict iV as a function of the initial sea state. By ob- 
taining N for each sea state, we can quantitatively predict 
the enhancement of the freak wave occurrence probability 
for that sea state. 

However, this is not a completed investigation. After fully 
understanding the quantitative consequences of the nonlin- 
ear effect, the ultimate goal of the project is to combine 
nonlinearity and deflection by random currents in a single 
model, which will allow the probability of rogue wave for- 
mation to be predicted for a wide range of realistic sea con- 
ditions. Previous linear results show that the effect of cur- 
rent is well characterized by a single parameter: the freak 
index [Heller et al., 2008; Ying et ai, 2011]. Preliminary 
results of combining nonlinearity and deflection by random 
currents show that the freak index is still a dominant pa- 
rameter. 

When we introduce a nonzero current field, additional 
length scales come in to play, including the characteristic 
eddy size (the scale on with the random current varies) and 
the typical distance required for the first focal point to ap- 
pear. It is still unclear how these parameters will interact 
with the wavelength of the incoming wave and the typical 
distance scale for the buildup of the nonlinear effect. An 
in-depth investigation is required to understand the under- 
lying mechanism through which the formation of hot and 
cold spots is aided by nonlinear focusing. 

Also, there is a clear need to compare the model simula- 
tions with observations and experiments. Although compre- 
hensive global data are not available at this point, it may 
be possible to compare the results with local observations 
where data are more readily available, e.g. in the North 
Sea. 
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